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INTRODUCTION 

A  "simulation”  is  an  experiment  run  inside  the  compu¬ 
ter  (instead  of  in  what  is  generally,  in  the  field,  termed 
the  "real  world”)  in  order  to  obtain  information  about  a 
system  (existing  or  proposed),  modifications  to  a  system, 
or  several  competing  systems.  Since  it  is  well-known  that 
much  more  information  can  be  extracted  from  experiments 
which  are  carefully  designed  statistically  than  can  be 
from  experiments  which  are  not  statistically  designed, 
statistics  has  a  large  role  to  play  in  simulation.  This 
chapter  deals  with  that  role.  This  is  an  especially 
important  chapter  for  those  not  used  to  thinking  of  their 
simulations  in  statistical  terms  (e.g.  ,  those  who  "run 
until  the  money  runs  out"  or  who  believe  that  "one  run 
will  show  me  how  the  system  behaves"),  as  it  will  enable 


This  research  was  supported  by  Office  of  Naval  Research 
Contract  No.  N00014-78-C-0543. 


I 


«> 

<c 


than  for  valid  results  obtained  without  a  carefully 
designed  experiment ,  and  to  obtain  valid  results  for 
systems  where  (without  using  statistical  design)  the  study 
is  computationally  infeasible. 

As  we  have  said,  since  simulation  data  are  essenti¬ 
ally  obtained  from  an  experiment  run  "inside  the  computer," 
we  need  to  know  how  to  design  the  experiment  whereby  those 
data  are  generated,  and  we  need  to  know  how  to  analyze  the 
resulting  data.  Now  this  comprises  a  large  part  of  the 
field  commonly  termed  statist i cs ,  to  which  whole  university 
departments  are  devoted.  Hence,  only  the  most  important 
aspects  of  design  and  analysis  will  be  covered  here,  with 
a  cavaet  that,  if  a  simulation  study  is  to  be  properly 
designed  and  executed,  in  many  cases  it  will  be  desirable 
to  have  associated  with  the  study  from  the  outset  a 
statistician  competent  in  simulation  applications. 

It  is  a  fact  that  statistics  grew  up  in  a  largely 
agricultural  setting  where  such  assumptions  as  homoscedas- 
ticity  (equality  of  variances  of  observations  taken  from 
diverse  sources,  such  as  weights  of  tomotoes  from  differ¬ 
ent  types  of  tomato  plants)  commonly  held  true,  while  in 
simulation  these  are  often  violated  (e.g.,  serious  hetero- 
scadastieity  —  greatly  unequal  variances  --  is  more  often 
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the  rule  in  simulation  studies,  as  when  a  number  of 
diverse  alternative  job -shop  scheduling  rules  are  being 
simulated).  The  theory  required  to  handle  these  new 

situations  has  only  recently  been  developed  (e.g.,  for 
heteroscedasticity  [ll  and  for  autocorre  lated  observations 
[2]).  Thus,  many  parts  of  the  topics  covered  here  will, 
for  some  time  to  come,  be  unavailable  in  almost  any  texts 
or  books  and  manuals  on  statistical  methods. 

1 .  Design 

1 . a .  Factorial,  fractional -factorial ,  and  one-at-a- 
t ime  experiments:  why  and  why /not,  efficiency . 

The  pitfalls  of  lack  of  statistical  design  in  a 
simulation  experiment  include:  invalid  inferences;  valid 
inferences  at  substantially  increased  costs;  and  inabili¬ 
ty  to  complete  the  study. 

For  example,  suppose  a  simulation  model  has  been 
built  of  output  Y(x^,  x^)  of  a  system  which  depends  on  two 
input  variables  (temperature  in  °F.)  and  (reaction 
time  in  minutes)  [each  reader  will  know  similar  examples 
in  his/her  own  field,  and  may  find  it  helpful  to  re-cast 
this  example  in  the  terms  applicable  to  one  such  example]. 
Suppose  that  it  is  standard  to  operate  at  x^  =  300, 

Xj  =  10,  and  the  simulation  is  to  evaluate  possible  gains 
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of  increasing  x1  (by  25°F.)  or  x ^  (by  2  minutes). 
Commonly  one  would  run  the  model  at  (x^ ,  x^)  equal  to 

(300,  ]0),  (325,  10),  (300,  12)  to  evaluate  the  effect 
of  the  proposed  changes  ( one-at-a-tlme  method).  Suppose 
one  does  so,  and  the  results  (each  based  on  one  run)  are 
as  in  Table  L.  Let  us  suppose  that  no  variability  is 


Table  1.  Results  of  one-at--a-t ime  experiment. 


( X-J  )  x 2  ) 

Y(x^,  Xp)  from  simulation  run 

(300,10) 

100.0 

(325,  10) 

102.5 

(300,  12) 

101.0 

present  (so,  untypically,  results  of  one  simulation  run 
are  typical),  that  past  experience  with  the  system  has 
typically  yielded  100.0  output  units  at  (x^,  x0 )  =  (  300  , 
10),  that  expensive  plant  changes  are  required  to  imple¬ 
ment  changes  of  this  magnitude  (so  we  must  simulate  before 
making  changes  in  the  plant),  and  that  to  be  economically 
feasible  we  must  obtain  at  least  a  5%  process  yield  in¬ 
crease  from  any  recommended  changes.  We  estimate  the 
gain  from  incrementing  xq  to  be 

102.5  -  100.0  =  2.5 

and  the  gain  from  incrementing  x,  to  be 

L. 
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for  a  total  estimated  gain  of  2.5  +  1.0  =  3.5  units,  less 
than  5%,  so  the  proposed  improvement  in  process  is  dis¬ 
carded  . 

Now  the  above  assumes  an  output  linear  in  each  of  x^ 
and  (perhaps  reasonable  in  a  small  range  about  the 
usual  process  operating  conditions),  but  with  no  synergis¬ 
tic  effects  (no  interaction ) . . .  which  is  often  false  in 
the  real-world  (at  least,  we  should  not  believe  such  an 
assertion  without  validation,  and  the  design  above  allows 
no  such  verification).  For  example,  an  experiment  run  at 
(325,  12)  might  complete  our  data  set  to  that  of  Table  2. 
This  set  of  data  has  a  far  different  interpretation  than 


2 

Table  2.  Results  of  2  experiment. 


(x^,  X2) 

Y(x^,  x2 )  from  simulation  run 

(300,  10) 

100. 0 

(325,  10) 

102.  5 

(300,  12) 

101.0 

(325,  12) 

110.0 

_ -  „  ■  „  .  -  - —  —  ■  - - - - - -  ■  ■ 

does  that  of  Table  1;  we  now  see  it  will  be  very  profi¬ 
table  (double  the  needed  5%  gain  in  output)  to  make  the 
plant  modifications  needed  to  run  at  the  higher  temperature 
(x^  =  325°F)  and  longer  time  (X2  =  12  min.),  thus  avoiding 
an  invalid  inference.  The  results  of  Table  2  would  be 


2 

better  displayed  as  in  Table  3,  which  shows  the  2  (two 
factors,  each  at  two  levels)  design  more  clearly.  Such 
results  as  these  can  be  obtained  from  as  simple  an  under¬ 


lying  true  function  as 

Y(x1,x?>  =  100.0+0.1(x1-300.0)+0.5(x2-10.0) 
+0.2<x1-300.0)(xo-10.0). 

As  the  number  n  of  variables  increases,  the  (invalid) 
one-at-a-time  method  will  require  a  number  of  simulation 
runs  equal  to  the  number  of  variables;  e.g.  if  in 
addition  to  xx  (temperature)  and  x2  (reaction  time)  we 
also  have  variables  x^  (arid  concentration)  and  x8j  (pH) 
present,  then  n  =  h  simulation  runs  would  be  needed  by 
this  (invalid)  method.  In  order  to  assess  effects  with 

each  variable  at  the  traditional  or  a  higher  level,  the 

4 

analog  cf  the  experiment  in  Table  2  would  require  2-16 
simulation  runs.  As  these  runs  can  be  very  expensive 
(e.g.,  in  nuclear  simulations,  where  a  run  can  consume 
several  hours  of  computer  time),  we  now  wish  to 
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comparably  correct  inferences  at  substantially  reduced 

cost  ■ 

To  perform  this  task,  we  must  first  consider  models 
for  yield  Y  as  a  function  of  variables  x^  ,  x4  ,  .... 

It  is  reasonable  to  assume  that  true  mean  yield  E(Y)  may 
be  accurately  represented  by  a  polynomial  equation  of 
sufficiently  high  order: 


EY  ( x  2  , 

X2  ’  x3 ’ 

x4  > 

• 

.  .  ) 

=  60  + 

Blxl  + 

$2x2 

+ 

03x3  + 

^4X4 

+ 

612X1X2 

+  6 

13 

xlx  3  + 

614xlx4+s2  3x2x3+3.?4x:x4+e34x3x« 

+ 

811X1  + 

B22 

+  63  3x 

•? 

3  + 

644X4 

+ 

8123x1x 

2x3 

♦ 

8124x1x 

2X4 

+  8134x1x3x[4  +  8234X2X3X4 

+ 

81234x1 

x2x3 

X4 

+  ... 

* 

Different  choices  of  the  points  (x-p  x?,  x^,  x4 ,  ...)  at 
which  simulation  runs  should  be  made  (each  such  choice  is 
called  an  experimental  design)  allow  us  to  estimate  var¬ 
ious  of  the  g’s  (singly  or  in  combination),  while  assuming 
others  are  negligible  (which  can  often  be  tested). 
Typically  one  assumes  terms  of  higher  than  second  order 
are  negligible  (0  =  8123  =  8124  =  8134  =  8234  =  81234  = 
...).  The  coefficients  g.  ,  g,,,  03,  64  are  called  main 
ef  fects  of  variables  Xp  x2  ,  x3,  x4  ,  while  gp,  ®13’  ®14’ 
833,  ®24’  ^34  are  caHed  2- factor  interactions  (2  f.i.) 


of  the  variables  involved  in  the  subscripts.  An  experi¬ 
mental  design  is  called  a  Resolution  III  design  if  no 

main  effect  is  conf ounded  with  (not  able  to  be  estimated 
sepai ate]y  : tom)  any  rthei  main  effect,  but  main 
effects  are  confounded  with  2  f.i.  and  2  f.i.  with  each 
other.  In  a  Resolution  IV  design  no  main  effect  is  con¬ 
founded  with  any  other  main  effect  or  2  f.i. ,  but  2  f.i.'s 
are  confounded  with  each  other.  Thus,  Resolution  III  will 
allow  us  to  fit  a  model  Bq  ♦  l  ^ixj_  in  factors  ,  but 
interaction  will  bias  the  fit.  A  resolution  IV  design 
will  allow  a  fit  unbiased  by  2  f.i. ,  i.e.  will  estimate 
(not  +_  some  of  the  B^j's  of  interactions). 

The  one-at-a-time  method  yields  a  resolution  IV  design, 
and  its  disadvantages  have  been  noted  above. 

l.b.  Screening  experiments. 

In  screening  experiments  one  attempts  to  use  designs 
which  allow  one  to  find  out  (at  relatively  modest  cost, 
i.e.  no.  of  simulation  runs)  which  of  the  variables  , 
x2 ,  ,  x4,  ...  are  most  influential  as  to  one's  system's 

output  in  very  few  experiments  by  using  a  Resolution  III 
design  (of  course  one  then  gives  up  almost  all  model 
fitting  ability,  the  intention  being  to  run  a  more  exten¬ 
sive  experiment  later  on  the  few  highly  influential 
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variables  identified).  For  example,  one  design  commonly 
used  is  the  Resolution  III  Placket t-Burman  design  which 
can  study  7  variables  x, ,  x2  »  x3,  x4 ,  x&,  xfi,  Xj  in  B 

simulation  runs,  shovin  in  Table  u .  In  this  design. 


Table  4.  Resolution  III  Plackett -Burman  design. 


_ 

*1 

X2 

Variable 

x3  x4 

X5 

x6 

X? 

1 

4 

4 

4 

- 

4 

- 

- 

2 

- 

4 

4 

+ 

- 

4 

- 

Experiment 

3 

- 

- 

4 

4 

4 

- 

4 

4 

4 

- 

- 

4- 

4 

4 

- 

5 

- 

+ 

- 

- 

4 

+ 

+ 

6 

4 

- 

■f 

- 

- 

4 

4 

7 

8 

4 

+ 

4 

— 

4 

each  experimental  variable  has  2  levels,  low(-)  and  high(+) 
(e.g.,  3C0(->  and  325(+)  for  x^).  In  experiment  1,  varia¬ 
bles  1,  2,  3,  5  are  set  to  their  high  levels,  while  varia¬ 
bles  4,  G,  7  are  set  to  their  low  levels.  If  desired 
later,  this  design  can  be  augmented  to  a  Resolution  IV 
design  by  adding  8  more  runs  (for  a  total  of  16  runs).  One 
can  also  study  6  factors  in  a  12  experiment  (Webb)  experi¬ 
mental  design,  at  Resolution  IV.  As  another  example  of  a 
screening  design,  if  we  have  8  independent  variables  then 


a  Webb  design  of  Resolution  IV  is  available  (see  Table  5) 


which  involves  10  experiments. 


Table  5.  5-variable  screening  design. 


Variable 


x2  x5  x4 


Experiment 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 


+ 

+ 

+ 

+ 


+ 

+ 

+ 

+ 


+ 

+ 

+ 

+ 


Note  that  the  first  5  experiments  form  a  Resolution  III 
design,  but  the  10-run  experiment  is  much  more  efficient 
than  the  5-run  experiment  (it  yields  estimates  with 
variance  o2/9  for  main  effects,  vs_.  o2  for  the  5-run  one- 
at-a-time  experiment...  a  9-fold  reduction  at  a  price  of 
5  more  runs  ,  and  with  2  f.i.  elimination;  here  oz  is  the 
variability  inherent  in  each  simulation  run's  outcome). 

As  a  final  example  of  a  screening  design,  with  4  variables 
,  x2»  x3,  x4  a  Webb  Resolution  IV  design  with  8 
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experiments  is  available  (Table  6)  (which  is  a 


"fractional  factorial"  design,  i.e.  a  fraction  of  the 
2^  =  16  experiment  factorial  design  with  confounding 
scheme  I  =  ABCD)  which  allows  one  to  fit  Bg ,  >  63’ 

unbiased  by  two-factor  interactions  (though  those  inter 
actions  themselves  are  not  estimable). 


Table  6.  5-variable  screening  design. 


Variable 

*1  *2  x3  x4 

1 

2 

+  -  -  + 

3 

+  -  + 

4 

+  + 

Experiment  5 

+  + 

6 

+  -  + 

7 

+  4- 

8 

+  +  +  + 

1 . c .  Central  composite  designs  and  full  quadratic 
models . 

Assuming  screening  has  been  completed  we  will  wish 
to  use  a  design  which  allows  for  assessment  of  all  main 
ef f ects ,  all  2  f.i,,  and  (perhaps)  all  quadratic  effect s . 


We  will  illustrate  with  a  4  variable  example. 

Here  one  wishes  not  simply  to  find  the  most  important 
factors,  but  rather  to  model  system  output  as  a  function 
of  design  settings.  Of  course  the  appropriate  design  to 
use  depends  on  one’s  goals  and  scope  of  study,  as  well  as 
on  one's  budget.  One  possibility  is  an  11  experiment  Webb 
design  (Table  7).  This  design  allows  one  to  fit  , 

$2  >  $  3  >  ^4’  ^12  *  ^13’  ^14'  ^2  3 ’  ^24'  ^34  ® umin g  no 

curvature . 


Table  7  .  4-variable  Webb  design.* 


Experiment 


X1 

Variable 

x2  x3 

x4 

0 

0 

0 

0 

0 

0 

1 

1 

0 

1 

0 

1 

1 

0 

1 

1 

0 

1 

0 

0 

1 

0 

1 

0 

1 

0 

0 

*  "0"  denotes  low  (-),  "1"  denotes  high  (+) 
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Since  curvature  cannot  often  be  ruled  out  a  priori, 
one  will  usually  desire  another  design  (which,  however, 
will  require  more  experiments)  unless  one's  budget  is 
severely  constrained.  The  full  factorial  24  design  (16 
experiments )  allows  one  to  estimate  Bq  +  +  B22  +  B33  + 

B44’  Si»  ’  83  ’  84’  812  ’  813  *  814’  g23’  824’  834’  8123’  8124’ 
8134’  8 2 3 4  and  Assuming  B1234  negligible)  have  an  estimate 
of  experimental  error  a2.  [No  suitable  "fractions"  of  this 
design  exist,  as  they  all  confound  one-  and/or  two-factor 
effects,  so  no  such  easy  reduction  of  number  of  experiments 
needed  is  possible  in  the  4-factor  case. ]  By  adding  one 
experiment,  one  obtains  the  24+I  design  (17  points)  ,  with 
which  one  can  now  separately  estimate  SQ,  +  B29  +833  + 

,  hence  (barring  cancelling  magnitudes)  assess  the  total 
quadratic  effect  independently  of  the  response  at  the 
center  ( Bq ) • 

While  the  above  designs  are  in  common  use,  more 
recently  the  Central  Composite  Design  (CCD)  hoc  Men  used 

to  good  advantage  in  such  situations.  This  requires  8  star 

*  *  4 

(or  axial)  points  be  added  to  the  2  +1  design,  for  a  total 

of  25  experiments,  and  allows  a  full  quadratic  model  to  be 

fitted.  Suitable  fractions  of  the  24  may  be  used  (i.e.  a 


fraction  with  at  most  one  2-factor  interaction  in  any  alias 
set),  and  (via  confounding  scheme  I  =  ABC)  one  can  obtain 
the  full  quadratic  model  estimation  with 


^  +  1  ♦  8  *  17  points. 


As  an  example ,  suppose  we  have  six  variables  to  consider 
say  xL,  x2,  x3,  x4,  x5,  x6.  A  full-factorial  approach 
would  require  26  =  64  experiments,  which  in  many  cases 
would  not  be  feasible.  However,  a  CCD  can  be  implemented 
with  only  26"2  +  <2)<6)  +  1  =  29  experiments  as  shown  in 
Table  8 .  There  the  levels  are  "coded"  so  that  "1"  means 
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Table  8.  CCD  studying  6  variables,  in  29  experiments, 


Name 


Expt .  No. 


X1 

x2 

Variable 

*3  x4 

X5 

X6 

1 

-1 

-1 

1 

-1 

-1 

i 

—  1 

1 

-1 

-1 

-1 

-1 

1 

1 

-1 

-1 

1 

1 

1 

1 

-1 

-1 

1 

-1 

-I 

-1 

1 

-1 

-1 

1 

-1 

-1 

1 

-1 

-1 

-1 

1 

-1 

1 

-1 

1 

1 

1 

-1 

1 

-1 

1 

-1 

-1 

-1 

1 

-1 

1 

-1 

-1 

1 

-1 

-1 

1 

-1 

-1 

1 

1 

1 

1 

-1 

-1 

1 

Center 


the  "high"  level  of  the  variable,  ,r-l"  means  the  "low" 
level  of  the  variable,  and  "0"  means  the  average  of  the 
high  and  low  levels.  "-a"  and  "a"  represent  multiples  of 
the  "low"  and  "high"  levels;  e.g.,  if  one  takes  a  =  1 
(face-centered  star  points)  these  are  equal  to  the 
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respective  "low”  and  "high”  levels  of  the  variable  in 
question.  Some  typical  possibilities  are  giver,  in  Table 
5  . 

Table  1.  Coding  of  variable  x^  levels,  CCD. 


-a 

1 

Level 

0 

1 

a 

a  =  1 

300 

300 

312 . 5 

325 

325 

a  =  2 

287.5 

300 

312.5 

325 

337.5 

If  one  has  bounds  and  Ih  on  variable  x^  and  wishes  to 
explore  the  full  space ,  a  =  1  is  often  recommended.  In 
other  settings  1  <  a  <  2  is  often  used.  The  specifics 
vary  from  setting  to  setting,  with  a  =  1.5  being  a 
reasonable  compromise  for  experimenters  who  do  not  have 
access  to  a  statistical  design  expert. 

1 . d .  Efficient  implementation. 

It  is  important  to  note  that  the  ordering  of  experi¬ 
ments  in  the  above  tables  is  not  the  recommended  order  in 
which  the  simulation  runs  should  be  made.  Since  often  the 
"next"  run  starts  with  the  end  of  the  "last"  run's  random 
number  stream,  runs  should  be  in  random  order  to  prevent 

systematic  effects  (from  possible  deviations  from  random¬ 
ness  of  one's  random 
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number  generator)  from  systematically  biasing  the 
results.  A  carefully  tested  and  chosen  random  number 
generator  is  essential,  and  an  existing  generator  should 

not  be  used  without  extensive  validation  (available  in 

[3]). 

Note  that  (unless  one  has  prior  knowledge  of  how  his 
system  will  use  random  numbers  and  correlation,  which  is 
rare)  use  of  variance  reduction  techniques  (see  Chapter  D) 
will  not  often  be  appropriate. 

1 . e .  Two-stage  and  sequential  designs . 

The  designs  explored  so  far  are  reasonable  ones  to 
use  in  homoscedast ic  systems.  In  heteroscedast ic  systems 
(where  o2  is  a  function  of  the  levels  of  xlt  ,  ...), 
other  designs  are  called  for.  Most  commonly  one  finds 
that  if  *  x2 ’  represent  levels  of  a  continuous 

variable  (such  as  temperature),  then  the  designs  given 
thus  far  are  appropriate  even  if  o2  varies  slightly  as 
x.  ,  x,,  ,  ...  are  changed. 

However  if  x^  ,  ,  ...  represent  the  presence  (1)  or 

absence  (0)  of  an  attribute,  then  new  designs  for  select¬ 
ing  the  best  combination  of  attributes  re:  system 
performance  are  called  for.  See  [4]  for  some  details. 
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2 .  Analysis  . 

The  designs  presented  in  section  1  are  appropriate 
for  situations  where  one  wants  to  know  "which  of  k  factors 
are  most  important  in  determining  my  system's  output?"  or 
"What  model  for  system  output,  in  k  input  factors,  should 
be  subjected  to  optimization  (Chapter  F)?"  Design  and 
analysis  of  the  simulation  experiment  when  one's  goal  is 
not  to  answer  the  above  questions,  but  rather  others,  such 
as:  a.  What  is  the  long-run  mean  of  my  system?  (which 

may  be  answered  using  a  trans format  ion -based  analysis, 
and  involves  questions  of  normality  and  run-in  time,  and 
leads  us  to  two-stage  and  regenerative  approaches); 

b.  Which  system  parameters  have  significant  effects? 
(which  leads  us  into  analysis  of  variance  (ANOVA),  and 
where  transformations  should  generally  not  be  used); 

c .  How  different  are  the  various  systems  1  performances? 
(which  requires  simultaneous  interval  estimates);  and 

d .  Which  is  the  best  system  (or  set  of  system  parameters)? 
(which  requires  the  new  methodology  of  ranking-and- 
selection  procedures),  is  a  subject  now  undergoing  rapid 
development.  While  such  problems  have  traditionally  been 
approached  in  the  past  with  transformations,  one- stage 
procedures,  or  ANOVA,  the  pitfalls  of  some  of  these 
traditions  and  of  equal-sample-sizes,  and  the  new  methods 
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re  cent  ly  developed  which  should  be  used  in  the  luture  , 
are  explored  in  Cl3. 

Towards  the  future,  ongoing  developments  in  multi¬ 
variate  analogs  of  the  design  and  analysis  procedures 
presented  should  lead  to  future  procedures  which  are  arble 
to  simultaneously  consider  several  output  characteristics. 
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